ATOMFT: Solving ODES and DAEs Using Taylor Series

Taylor series methods compute a solution to an initial value problem in ordinary differential equations by expanding each component of the solution in a long Taylor series. The series terms are generated recursively using the techniques of automatic differentiation. The ATOMFT system includes a translator to transform statements of the system of ODES into a FORTRAN 77 object program that is compiled, linked with the ATOMFT runtime library, and run to solve the problem. We review the use of the ATOMFT system for nonstiff and stiff ODES, the propagation of global errors, and applications to differential algebraic equations arising from certain control problems, to boundary value problems, to numerical quadrature, and to delay problems. 1. TAYLOR SERIES METHOD The solution of an initial value problem in ordinary differential equations expanded as a Taylor series has been given as both a classical and a numerical method for many years. The work of Sir Isaac Newton contains a four term series expansion for a nonelementary ordinary differential equation. In 1946, Miller [l] used recurrence relations to compute Taylor series terms for the Airy integral. Others [2,3] have written translator programs using automatic differentiation to write object programs for solving ODES. Moore [4] solved ODES, evaluating the Taylor series remainder term in interval arithmetic to compute a guaranteed enclosure of the solution. Lohner [5] is the latest of many who have advanced Moore’s ideas for interval enclosures of Taylor series solutions for ODES. Rall [S] g ives other applications of Taylor series methods. The philosophy of the Taylor series method is totally different from that of other methods in the solution of ODES. We use a power series for the solution function that is very long compared to the usual fourth-order or twelfth-order methods. For an ODE whose solution is f(t), the series terms for f(t) expanded at the solution point with an arbitrary stepsize h and stored as reduced derivatives, F(n + 1) := F(n) 2. These reduced derivatives are the Taylor series terms. We calculate them up to the 30th term and beyond. With the long Taylor series, it is then possible to calculate the radius of convergence. This is the principal departure from other methods. The arbitrary stepsize h is adjusted to an optimum stepsize after the radius of convergence has been calculated. To properly control the local truncation error, the optimum stepsize is determined from the series length, the radius The authors would like to express their gratitude to R. Morris for the initial design and coding of the translator program, to J. Fauss, D. Lowery and M. Prieto for their work on series analysis, to R. Moore, M. Tabor, J. Weiss for many helpful suggestions, to R. Stanford, P. Breckheimer and K. Berryman at Jet Propulsion Labs for requesting user defined functions, and to J. Wright for the many bugs found.